MODULE ModAbund

!!$ This module contains the subroutine for calculating abundances as a function of line of sight column density to the star.
!!$ Programming conventions are from Nyhoff and Leestma, _Fortran 90 for Scientists and Engineers_, 1997
!!$
!!$ STH: Samuel Harrold

!!$ Declare modules to USE
  USE ModParms
  USE ModFuncs

!!$ Do not assume undeclared variables
  IMPLICIT NONE



CONTAINS



!!$ CalcAbund
!!$ Estimate abundances for OI, CII, e, H, H2
!!$ TODO: include chemistry
!!$ Source: Glassgold et al 2004, Figs 4, 5.
!!$ Note: Only valid for 10^19 <= ngas <= 10^24
!!$ Abundances are estimated at vertical column density in a disk. From communication with Lacy (6/21/11), column density to the star is about 10 times the vertical column density. Column density to the star is common parameter between pisco and other models in spite of different spacial dependencies of other physical parameters. Abundances are not forced to go to zero. Abundances are listed in order of importance. Abundance*number density of gas = number density of species.

  SUBROUTINE CalcAbund(Ntostar, AbundCII, AbundOI, Abunde, AbundH, AbundH2)

!!$ Declare variables
!!$ Specifiy input and output
!!$ Ntostar = column density of gas in line of sight to star in particles cm^-2
!!$ Abund(OI/CII/e/H/H2) = abundances of OI/CII/e/H/H2 relative to all gas particles determined by column density to the star
    REAL, INTENT(IN)  :: Ntostar
    REAL, INTENT(OUT), OPTIONAL :: AbundCII, AbundOI, Abunde, AbundH, AbundH2

!!$ Specify internal variables
!!$ NAarray(CII/OI/e/H/H2) = arrays of coordinate pairs of 
!!$ column density to star (N) and abundance (A)
!!$ ((N1, A1),(N2, A2),...(Nn, An)) for CII and all collisional species. Approximated from Glassgold et al 2004. Note: only valid for 10^19 <= ngas <= 10^24
    REAL, DIMENSION(2, 2) :: NAarrayCII
    REAL, DIMENSION(3, 2) :: NAarrayOI
    REAL, DIMENSION(4, 2) :: NAarraye
    REAL, DIMENSION(2, 2) :: NAarrayH
    REAL, DIMENSION(5, 2) :: NAarrayH2

!!$ InterpLog interpolates a logarithmic function given two coordinate pairs:
!!$ y/y0 = (x/x0) ^ ( Log(y1/y0) / Log(x1/x0) )
!!$ TODO: make NAarrays allocateable. Make function that generates Abund for given data.

!!$ Extra/interpolate abundance of CII
    IF (PRESENT(AbundCII)) THEN
       NAarrayCII = RESHAPE( (/ &
            10.*8.E+20, 3.E-4, &
            10.*5.E+21, 1.E-9 /), (/2, 2/), &
            ORDER = (/2, 1/) )
       IF (Ntostar < NAarrayCII(1, 1)) THEN
          AbundCII =  NAarrayCII(1, 2)
       ELSE IF (NAarrayCII(1, 1) <= Ntostar) THEN
          AbundCII = InterpLog(x = Ntostar, &
               xi = NAarrayCII(1, 1), yi = NAarrayCII(1, 2), &
               xf = NAarrayCII(2, 1), yf = NAarrayCII(2, 2))
       END IF
    END IF

!!$ Extra/interpolate abundance of OI
    IF (PRESENT(AbundOI)) THEN
       NAarrayOI = RESHAPE( (/ &
            10.*1.E+21, 6.E-4, &
            10.*1.E+23, 2.E-5, &
            10.*3.E+23, 2.E-11 /), (/3, 2/), &
            ORDER = (/2, 1/) )
       IF (Ntostar < NAarrayOI(1, 1)) THEN
          AbundOI =  NAarrayOI(1, 2)
       ELSE IF (NAarrayOI(1, 1) <= Ntostar .AND. Ntostar <= NAarrayOI(2, 1)) THEN
          AbundOI = InterpLog(x = Ntostar, &
               xi = NAarrayOI(1, 1), yi = NAarrayOI(1, 2), &
               xf = NAarrayOI(2, 1), yf = NAarrayOI(2, 2))
       ELSE IF (NAarrayOI(2, 1) <= Ntostar) THEN
          AbundOI = InterpLog(x = Ntostar, &
               xi = NAarrayOI(2, 1), yi = NAarrayOI(2, 2), &
               xf = NAarrayOI(3, 1), yf = NAarrayOI(3, 2))
       END IF
    END IF

!!$ Extra/interpolate abundance of e
    IF (PRESENT(Abunde)) THEN
       NAarraye = RESHAPE( (/ &
            10.*1.E+18, 3.E-2, &
            10.*9.E+20, 3.E-4, &
            10.*4.E+21, 8.E-7, &
            10.*1.E+24, 3.E-9 /), (/4, 2/), &
            ORDER = (/2, 1/) )
       IF (Ntostar < NAarraye(1, 1)) THEN
          Abunde =  NAarraye(1, 2)
       ELSE IF (NAarraye(1, 1) <= Ntostar .AND. Ntostar < NAarraye(2, 1)) THEN
          Abunde = InterpLog(x = Ntostar, &
               xi = NAarraye(1, 1), yi = NAarraye(1, 2), &
               xf = NAarraye(2, 1), yf = NAarraye(2, 2))
       ELSE IF (NAarraye(2, 1) <= Ntostar .AND. Ntostar < NAarraye(3, 1)) THEN
          Abunde = InterpLog(x = Ntostar, &
               xi = NAarraye(2, 1), yi = NAarraye(2, 2), &
               xf = NAarraye(3, 1), yf = NAarraye(3, 2))
       ELSE IF (NAarraye(3, 1) <= Ntostar) THEN
          Abunde = InterpLog(x = Ntostar, &
               xi = NAarraye(3, 1), yi = NAarraye(3, 2), &
               xf = NAarraye(4, 1), yf = NAarraye(4, 2))
       END IF
    END IF

!!$ Extra/interpolate abundance of H
    IF (PRESENT(AbundH)) THEN
       NAarrayH = RESHAPE( (/ &
            10.*4.E+21, 1., &
            10.*1.E+24, 1.E-4 /), (/2, 2/), &
            ORDER = (/2, 1/) )
       IF (Ntostar < NAarrayH(1, 1)) THEN
          AbundH =  NAarrayH(1, 2)
       ELSE IF (NAarrayH(1, 1) <= Ntostar) THEN
          AbundH = InterpLog(x = Ntostar, &
               xi = NAarrayH(1, 1), yi = NAarrayH(1, 2), &
               xf = NAarrayH(2, 1), yf = NAarrayH(2, 2))
       END IF
    END IF

!!$ Extra/interpolate abundance of H2
    IF (PRESENT(AbundH2)) THEN
       NAarrayH2 = RESHAPE( (/ &
            10.*1.0E+18, 3.E-5, &
            10.*7.0E+20, 1.E-6, &
            10.*1.2E+21, 2.E-3, &
            10.*3.0E+21, 1.E-3, &
            10.*7.0E+21, 5.E-1 /), (/5, 2/), &
            ORDER = (/2, 1/) )
       IF (Ntostar < NAarrayH2(1, 1)) THEN
          AbundH2 =  NAarrayH2(1, 2)
       ELSE IF (NAarrayH2(1, 1) <= Ntostar .AND. Ntostar < NAarrayH2(2, 1)) THEN
          AbundH2 = InterpLog(x = Ntostar, &
               xi = NAarrayH2(1, 1), yi = NAarrayH2(1, 2), &
               xf = NAarrayH2(2, 1), yf = NAarrayH2(2, 2))
       ELSE IF (NAarrayH2(2, 1) <= Ntostar .AND. Ntostar < NAarrayH2(3, 1)) THEN
          AbundH2 = InterpLog(x = Ntostar, &
               xi = NAarrayH2(2, 1), yi = NAarrayH2(2, 2), &
               xf = NAarrayH2(3, 1), yf = NAarrayH2(3, 2))
       ELSE IF (NAarrayH2(3, 1) <= Ntostar .AND. Ntostar < NAarrayH2(4, 1)) THEN
          AbundH2 = InterpLog(x = Ntostar, &
               xi = NAarrayH2(3, 1), yi = NAarrayH2(3, 2), &
               xf = NAarrayH2(4, 1), yf = NAarrayH2(4, 2))
       ELSE IF (NAarrayH2(4, 1) <= Ntostar .AND. Ntostar < NAarrayH2(5, 1)) THEN
          AbundH2 = InterpLog(x = Ntostar, &
               xi = NAarrayH2(4, 1), yi = NAarrayH2(4, 2), &
               xf = NAarrayH2(5, 1), yf = NAarrayH2(5, 2))
       ELSE IF (NAarrayH2(5, 1) <= Ntostar) THEN
          AbundH2 = NAarrayH2(5, 2)
       END IF
    END IF

  END SUBROUTINE CalcAbund



END MODULE ModAbund
